Tumor-immune metaphenotypes orchestrate an evolutionary bottleneck that promotes metabolic transformation.

Introduction: Metabolism plays a complex role in the evolution of cancerous tumors, including inducing a multifaceted effect on the immune system to aid immune escape. Immune escape is, by definition, a collective phenomenon by requiring the presence of two cell types interacting in close proximity: tumor and immune. The microenvironmental context of these interactions is influenced by the dynamic process of blood vessel growth and remodelling, creating heterogeneous patches of well-vascularized tumor or acidic niches. Methods: Here, we present a multiscale mathematical model that captures the phenotypic, vascular, microenvironmental, and spatial heterogeneity which shapes acid-mediated invasion and immune escape over a biologically-realistic time scale. The model explores several immune escape mechanisms such as i) acid inactivation of immune cells, ii) competition for glucose, and iii) inhibitory immune checkpoint receptor expression (PD-L1). We also explore the efficacy of anti-PD-L1 and sodium bicarbonate buffer agents for treatment. To aid in understanding immune escape as a collective cellular phenomenon, we define immune escape in the context of six collective phenotypes (termed "meta-phenotypes"): Self-Acidify, Mooch Acid, PD-L1 Attack, Mooch PD-L1, Proliferate Fast, and Starve Glucose. Results: Fomenting a stronger immune response leads to initial benefits (additional cytotoxicity), but this advantage is offset by increased cell turnover that leads to accelerated evolution and the emergence of aggressive phenotypes. This creates a bimodal therapy landscape: either the immune system should be maximized for complete cure, or kept in check to avoid rapid evolution of invasive cells. These constraints are dependent on heterogeneity in vascular context, microenvironmental acidification, and the strength of immune response. Discussion: This model helps to untangle the key constraints on evolutionary costs and benefits of three key phenotypic axes on tumor invasion and treatment: acid-resistance, glycolysis, and PD-L1 expression. The benefits of concomitant anti-PD-L1 and buffer treatments is a promising treatment strategy to limit the adverse effects of immune escape.

Modeling phenotypic heterogeneity towards evolutionarily inspired osteosarcoma therapy.

Osteosarcoma is the most common bone sarcoma in children and young adults. While universally delivered, chemotherapy only benefits roughly half of patients with localized disease. Increasingly, intratumoral heterogeneity is recognized as a source of therapeutic resistance. In this study, we develop and evaluate an in vitro model of osteosarcoma heterogeneity based on phenotype and genotype. Cancer cell populations vary in their environment-specific growth rates and in their sensitivity to chemotherapy. We present the genotypic and phenotypic characterization of an osteosarcoma cell line panel with a focus on co-cultures of the most phenotypically divergent cell lines, 143B and SAOS2. Modest environmental (pH, glutamine) or chemical perturbations dramatically shift the success and composition of cell lines. We demonstrate that in nutrient rich culture conditions 143B outcompetes SAOS2. But, under nutrient deprivation or conventional chemotherapy, SAOS2 growth can be favored in spheroids. Importantly, when the simplest heterogeneity state is evaluated, a two-cell line coculture, perturbations that affect the faster growing cell line have only a modest effect on final spheroid size. Thus the only evaluated therapies to eliminate the spheroids were by switching therapies from a first strike to a second strike. This extensively characterized, widely available system, can be modeled and scaled to allow for improved strategies to anticipate resistance in osteosarcoma due to heterogeneity.

Virtual alignment of pathology image series for multi-gigapixel whole slide images.

Interest in spatial omics is on the rise, but generation of highly multiplexed images remains challenging, due to cost, expertise, methodical constraints, and access to technology. An alternative approach is to register collections of whole slide images (WSI), generating spatially aligned datasets. WSI registration is a two-part problem, the first being the alignment itself and the second the application of transformations to huge multi-gigapixel images. To address both challenges, we developed Virtual Alignment of pathoLogy Image Series (VALIS), software which enables generation of highly multiplexed images by aligning any number of brightfield and/or immunofluorescent WSI, the results of which can be saved in the ome.tiff format. Benchmarking using publicly available datasets indicates VALIS provides state-of-the-art accuracy in WSI registration and 3D reconstruction. Leveraging existing open-source software tools, VALIS is written in Python, providing a free, fast, scalable, robust, and easy-to-use pipeline for registering multi-gigapixel WSI, facilitating downstream spatial analyses.

Intermetastatic and Intrametastatic Heterogeneity Shapes Adaptive Therapy Cycling Dynamics.

Adaptive therapies that alternate between drug applications and drug-free vacations can exploit competition between sensitive and resistant cells to maximize the time to progression. However, optimal dosing schedules depend on the properties of metastases, which are often not directly measurable in clinical practice. Here, we proposed a framework for estimating features of metastases through tumor response dynamics during the first adaptive therapy treatment cycle. Longitudinal prostate-specific antigen (PSA) levels in 16 patients with metastatic castration-resistant prostate cancer undergoing adaptive androgen deprivation treatment were analyzed to investigate relationships between cycle dynamics and clinical variables such as Gleason score, the change in the number of metastases over a cycle, and the total number of cycles over the course of treatment. The first cycle of adaptive therapy, which consists of a response period (applying therapy until 50% PSA reduction), and a regrowth period (removing treatment until reaching initial PSA levels), delineated several features of the computational metastatic system: larger metastases had longer cycles; a higher proportion of drug-resistant cells slowed the cycles; and a faster cell turnover rate sped up drug response time and slowed regrowth time. The number of metastases did not affect cycle times, as response dynamics were dominated by the largest tumors rather than the aggregate. In addition, systems with higher intermetastasis heterogeneity responded better to continuous therapy and correlated with dynamics from patients with high or low Gleason scores. Conversely, systems with higher intrametastasis heterogeneity responded better to adaptive therapy and correlated with dynamics from patients with intermediate Gleason scores. SIGNIFICANCE: Multiscale mathematical modeling combined with biomarker dynamics during adaptive therapy helps identify underlying features of metastatic cancer to inform treatment decisions.

Adaptive therapy for ovarian cancer: An integrated approach to PARP inhibitor scheduling.

Toxicity and emerging drug resistance are important challenges in PARP inhibitor (PARPi) treatment of ovarian cancer. Recent research has shown that evolutionary-inspired treatment algorithms which adapt treatment to the tumor's treatment response (adaptive therapy) can help to mitigate both. Here, we present a first step in developing an adaptive therapy protocol for PARPi treatment by combining mathematical modelling and wet-lab experiments to characterize the cell population dynamics under different PARPi schedules. Using data from in vitro Incucyte Zoom time-lapse microscopy experiments and a step-wise model selection process we derive a calibrated and validated ordinary differential equation model, which we then use to test different plausible adaptive treatment schedules. Our model can accurately predict the in vitro treatment dynamics, even to new schedules, and suggests that treatment modifications need to be carefully timed, or one risks losing control over tumour growth, even in the absence of any resistance. This is because our model predicts that multiple rounds of cell division are required for cells to acquire sufficient DNA damage to induce apoptosis. As a result, adaptive therapy algorithms that modulate treatment but never completely withdraw it are predicted to perform better in this setting than strategies based on treatment interruptions. Pilot experiments in vivo confirm this conclusion. Overall, this study contributes to a better understanding of the impact of scheduling on treatment outcome for PARPis and showcases some of the challenges involved in developing adaptive therapies for new treatment settings.

Agent-based methods facilitate integrative science in cancer.

In this opinion, we highlight agent-based modeling as a key tool for exploration of cell-cell and cell-environment interactions that drive cancer progression, therapeutic resistance, and metastasis. These biological phenomena are particularly suited to be captured at the cell-scale resolution possible only within agent-based or individual-based mathematical models. These modeling approaches complement experimental work (in vitro and in vivo systems) through parameterization and data extrapolation but also feed forward to drive new experiments that test model-generated predictions.

Mistic: An open-source multiplexed image t-SNE viewer.

Understanding the complex ecology of a tumor tissue and the spatiotemporal relationships between its cellular and microenvironment components is becoming a key component of translational research, especially in immuno-oncology. The generation and analysis of multiplexed images from patient samples is of paramount importance to facilitate this understanding. Here, we present Mistic, an open-source multiplexed image t-SNE viewer that enables the simultaneous viewing of multiple 2D images rendered using multiple layout options to provide an overall visual preview of the entire dataset. In particular, the positions of the images can be t-SNE or UMAP coordinates. This grouped view of all images allows an exploratory understanding of the specific expression pattern of a given biomarker or collection of biomarkers across all images, helps to identify images expressing a particular phenotype, and can help select images for subsequent downstream analysis. Currently, there is no freely available tool to generate such image t-SNEs.

Spatial structure impacts adaptive therapy by shaping intra-tumoral competition.

Background: Adaptive therapy aims to tackle cancer drug resistance by leveraging resource competition between drug-sensitive and resistant cells. Here, we present a theoretical study of intra-tumoral competition during adaptive therapy, to investigate under which circumstances it will be superior to aggressive treatment. Methods: We develop and analyse a simple, 2-D, on-lattice, agent-based tumour model in which cells are classified as fully drug-sensitive or resistant. Subsequently, we compare this model to its corresponding non-spatial ordinary differential equation model, and fit it to longitudinal prostate-specific antigen data from 65 prostate cancer patients undergoing intermittent androgen deprivation therapy following biochemical recurrence. Results: Leveraging the individual-based nature of our model, we explicitly demonstrate competitive suppression of resistance during adaptive therapy, and examine how different factors, such as the initial resistance fraction or resistance costs, alter competition. This not only corroborates our theoretical understanding of adaptive therapy, but also reveals that competition of resistant cells with each other may play a more important role in adaptive therapy in solid tumours than was previously thought. To conclude, we present two case studies, which demonstrate the implications of our work for: (i) mathematical modelling of adaptive therapy, and (ii) the intra-tumoral dynamics in prostate cancer patients during intermittent androgen deprivation treatment, a precursor of adaptive therapy. Conclusion: Our work shows that the tumour's spatial architecture is an important factor in adaptive therapy and provides insights into how adaptive therapy leverages both inter- and intra-specific competition to control resistance.

Immunosuppressive niche engineering at the onset of human colorectal cancer.

The evolutionary dynamics of tumor initiation remain undetermined, and the interplay between neoplastic cells and the immune system is hypothesized to be critical in transformation. Colorectal cancer (CRC) presents a unique opportunity to study the transition to malignancy as pre-cancers (adenomas) and early-stage cancers are frequently resected. Here, we examine tumor-immune eco-evolutionary dynamics from pre-cancer to carcinoma using a computational model, ecological analysis of digital pathology data, and neoantigen prediction in 62 patient samples. Modeling predicted recruitment of immunosuppressive cells would be the most common driver of transformation. As predicted, ecological analysis reveals that progressed adenomas co-localized with immunosuppressive cells and cytokines, while benign adenomas co-localized with a mixed immune response. Carcinomas converge to a common immune "cold" ecology, relaxing selection against immunogenicity and high neoantigen burdens, with little evidence for PD-L1 overexpression driving tumor initiation. These findings suggest re-engineering the immunosuppressive niche may prove an effective immunotherapy in CRC.

Roadmap on plasticity and epigenetics in cancer.

The role of plasticity and epigenetics in shaping cancer evolution and response to therapy has taken center stage with recent technological advances including single cell sequencing. This roadmap article is focused on state-of-the-art mathematical and experimental approaches to interrogate plasticity in cancer, and addresses the following themes and questions: is there a formal overarching framework that encompasses both non-genetic plasticity and mutation-driven somatic evolution? How do we measure and model the role of the microenvironment in influencing/controlling non-genetic plasticity? How can we experimentally study non-genetic plasticity? Which mathematical techniques are required or best suited? What are the clinical and practical applications and implications of these concepts?

Fluctuating methylation clocks for cell lineage tracing at high temporal resolution in human tissues.

Molecular clocks that record cell ancestry mutate too slowly to measure the short-timescale dynamics of cell renewal in adult tissues. Here, we show that fluctuating DNA methylation marks can be used as clocks in cells where ongoing methylation and demethylation cause repeated 'flip-flops' between methylated and unmethylated states. We identify endogenous fluctuating CpG (fCpG) sites using standard methylation arrays and develop a mathematical model to quantitatively measure human adult stem cell dynamics from these data. Small intestinal crypts were inferred to contain slightly more stem cells than the colon, with slower stem cell replacement in the small intestine. Germline APC mutation increased the number of replacements per crypt. In blood, we measured rapid expansion of acute leukemia and slower growth of chronic disease. Thus, the patterns of human somatic cell birth and death are measurable with fluctuating methylation clocks (FMCs).

The role of memory in non-genetic inheritance and its impact on cancer treatment resistance.

Intra-tumour heterogeneity is a leading cause of treatment failure and disease progression in cancer. While genetic mutations have long been accepted as a primary mechanism of generating this heterogeneity, the role of phenotypic plasticity is becoming increasingly apparent as a driver of intra-tumour heterogeneity. Consequently, understanding the role of this plasticity in treatment resistance and failure is a key component of improving cancer therapy. We develop a mathematical model of stochastic phenotype switching that tracks the evolution of drug-sensitive and drug-tolerant subpopulations to clarify the role of phenotype switching on population growth rates and tumour persistence. By including cytotoxic therapy in the model, we show that, depending on the strategy of the drug-tolerant subpopulation, stochastic phenotype switching can lead to either transient or permanent drug resistance. We study the role of phenotypic heterogeneity in a drug-resistant, genetically homogeneous population of non-small cell lung cancer cells to derive a rational treatment schedule that drives population extinction and avoids competitive release of the drug-tolerant sub-population. This model-informed therapeutic schedule results in increased treatment efficacy when compared against periodic therapy, and, most importantly, sustained tumour decay without the development of resistance.

Normal tissue architecture determines the evolutionary course of cancer.

Cancer growth can be described as a caricature of the renewal process of the tissue of origin, where the tissue architecture has a strong influence on the evolutionary dynamics within the tumor. Using a classic, well-studied model of tumor evolution (a passenger-driver mutation model) we systematically alter spatial constraints and cell mixing rates to show how tissue structure influences functional (driver) mutations and genetic heterogeneity over time. This approach explores a key mechanism behind both inter-patient and intratumoral tumor heterogeneity: competition for space. Time-varying competition leads to an emergent transition from Darwinian premalignant growth to subsequent invasive neutral tumor growth. Initial spatial constraints determine the emergent mode of evolution (Darwinian to neutral) without a change in cell-specific mutation rate or fitness effects. Driver acquisition during the Darwinian precancerous stage may be modulated en route to neutral evolution by the combination of two factors: spatial constraints and limited cellular mixing. These two factors occur naturally in ductal carcinomas, where the branching topology of the ductal network dictates spatial constraints and mixing rates.

The harsh microenvironment in early breast cancer selects for a Warburg phenotype.

The harsh microenvironment of ductal carcinoma in situ (DCIS) exerts strong evolutionary selection pressures on cancer cells. We hypothesize that the poor metabolic conditions near the ductal center foment the emergence of a Warburg Effect (WE) phenotype, wherein cells rapidly ferment glucose to lactic acid, even in normoxia. To test this hypothesis, we subjected low-glycolytic breast cancer cells to different microenvironmental selection pressures using combinations of hypoxia, acidosis, low glucose, and starvation for many months and isolated single clones for metabolic and transcriptomic profiling. The two harshest conditions selected for constitutively expressed WE phenotypes. RNA sequencing analysis of WE clones identified the transcription factor KLF4 as potential inducer of the WE phenotype. In stained DCIS samples, KLF4 expression was enriched in the area with the harshest microenvironmental conditions. We simulated in vivo DCIS phenotypic evolution using a mathematical model calibrated from the in vitro results. The WE phenotype emerged in the poor metabolic conditions near the necrotic core. We propose that harsh microenvironments within DCIS select for a WE phenotype through constitutive transcriptional reprogramming, thus conferring a survival advantage and facilitating further growth and invasion.

Turnover Modulates the Need for a Cost of Resistance in Adaptive Therapy.

Adaptive therapy seeks to exploit intratumoral competition to avoid, or at least delay, the emergence of therapy resistance in cancer. Motivated by promising results in prostate cancer, there is growing interest in extending this approach to other neoplasms. As such, it is urgent to understand the characteristics of a cancer that determine whether or not it will respond well to adaptive therapy. A plausible candidate for such a selection criterion is the fitness cost of resistance. In this article, we study a general, but simple, mathematical model to investigate whether the presence of a cost is necessary for adaptive therapy to extend the time to progression beyond that of a standard-of-care continuous therapy. Tumor cells were divided into sensitive and resistant populations and we model their competition using a system of two ordinary differential equations based on the Lotka-Volterra model. For tumors close to their environmental carrying capacity, a cost was not required. However, for tumors growing far from carrying capacity, a cost may be required to see meaningful gains. Notably, it is important to consider cell turnover in the tumor, and we discuss its role in modulating the impact of a resistance cost. To conclude, we present evidence for the predicted cost-turnover interplay in data from 67 patients with prostate cancer undergoing intermittent androgen deprivation therapy. Our work helps to clarify under which circumstances adaptive therapy may be beneficial and suggests that turnover may play an unexpectedly important role in the decision-making process. SIGNIFICANCE: Tumor cell turnover modulates the speed of selection against drug resistance by amplifying the effects of competition and resistance costs; as such, turnover is an important factor in resistance management via adaptive therapy.See related commentary by Strobl et al., p. 811.

Searching for Goldilocks: How Evolution and Ecology Can Help Uncover More Effective Patient-Specific Chemotherapies.

Deaths from cancer are mostly due to metastatic disease that becomes resistant to therapy. A mainstay treatment for many cancers is chemotherapy, for which the dosing strategy is primarily limited by patient toxicity. While this MTD approach builds upon the intuitively appealing principle that maximum therapeutic benefit is achieved by killing the largest possible number of cancer cells, there is increasing evidence that moderation might allow host-specific features to contribute to success. We believe that a "Goldilocks Window" of submaximal chemotherapy will yield improved overall outcomes. This window combines the complex interplay of cancer cell death, immune activity, emergence of chemoresistance, and metastatic dissemination. These multiple activities driven by chemotherapy have tradeoffs that depend on the specific agents used as well as their dosing levels and schedule. Here we present evidence supporting the idea that MTD may not always be the best approach and offer suggestions toward a more personalized treatment regime that integrates insights into patient-specific eco-evolutionary dynamics.

Evolutionary dynamics of neoantigens in growing tumors.

Cancers accumulate mutations that lead to neoantigens, novel peptides that elicit an immune response, and consequently undergo evolutionary selection. Here we establish how negative selection shapes the clonality of neoantigens in a growing cancer by constructing a mathematical model of neoantigen evolution. The model predicts that, without immune escape, tumor neoantigens are either clonal or at low frequency; hypermutated tumors can only establish after the evolution of immune escape. Moreover, the site frequency spectrum of somatic variants under negative selection appears more neutral as the strength of negative selection increases, which is consistent with classical neutral theory. These predictions are corroborated by the analysis of neoantigen frequencies and immune escape in exome and RNA sequencing data from 879 colon, stomach and endometrial cancers.

Hybrid Automata Library: A flexible platform for hybrid modeling with real-time visualization.

The Hybrid Automata Library (HAL) is a Java Library developed for use in mathematical oncology modeling. It is made of simple, efficient, generic components that can be used to model complex spatial systems. HAL's components can broadly be classified into: on- and off-lattice agent containers, finite difference diffusion fields, a GUI building system, and additional tools and utilities for computation and data collection. These components are designed to operate independently and are standardized to make them easy to interface with one another. As a demonstration of how modeling can be simplified using our approach, we have included a complete example of a hybrid model (a spatial model with interacting agent-based and PDE components). HAL is a useful asset for researchers who wish to build performant 1D, 2D and 3D hybrid models in Java, while not starting entirely from scratch. It is available on GitHub at https://github.com/MathOnco/HAL under the MIT License. HAL requires the Java JDK version 1.8 or later to compile and run the source code.

Correction: Acidity promotes tumour progression by altering macrophage phenotype in prostate cancer.

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Model genotype-phenotype mappings and the algorithmic structure of evolution.

Cancers are complex dynamic systems that undergo evolution and selection. Personalized medicine approaches in the clinic increasingly rely on predictions of tumour response to one or more therapies; these predictions are complicated by the inevitable evolution of the tumour. Despite enormous amounts of data on the mutational status of cancers and numerous therapies developed in recent decades to target these mutations, many of these treatments fail after a time due to the development of resistance in the tumour. The emergence of these resistant phenotypes is not easily predicted from genomic data, since the relationship between genotypes and phenotypes, termed the genotype-phenotype (GP) mapping, is neither injective nor functional. We present a review of models of this mapping within a generalized evolutionary framework that takes into account the relation between genotype, phenotype, environment and fitness. Different modelling approaches are described and compared, and many evolutionary results are shown to be conserved across studies despite using different underlying model systems. In addition, several areas for future work that remain understudied are identified, including plasticity and bet-hedging. The GP-mapping provides a pathway for understanding the potential routes of evolution taken by cancers, which will be necessary knowledge for improving personalized therapies.